GPS Footprinting

2020-05-26 19:51:48

Author: Cedric Huchuan Xia (email, github)

Affiliation: Penn Lifespan Informatics and Neuroimaging Center (PennLINC)

---

This GPS Footprinting project is inspired by Finn et al. Nat Neuro (2015) and Kaufmann et al. Nat Neuro (2015). These authors studied how personal differences in functional brain connectivity can identify individuals, mature during development, alter in neuropsychiatric illness, and differ between genders. The authors referred to these individual differences as brain fingerprinting.

Here, we apply the fingerprinting technique to highly sampled GPS data in a clinical sample of youth. We are interested in how their individual mobility patterns can distinguish one another’s identity, differ between genders, and alter in psychopathological groups. In other words, can our footprint tell us apart and something about our behavior?

GPS data preprocessing was performed according to Ian Barnett’s imputation algorithm, originally published in Biostatistics (2020).

1. Setup Environment

require(ggplot2)
require(summarytools)
require(cowplot)
require(caret)
require(corrplot)
require(RColorBrewer)
require(vembedr)
require(Rmisc)
require(varian)
source('~/Documents/GitHub/SmartphoneSensorPipeline/Extra/plotting_functions.R')

project_path = "~/Documents/xia_gps/"
data_path = file.path(project_path,"beiwe_output_043020")

gps_df_path = file.path(data_path,"Processed_Data/Group/feature_matrix.txt")

2. GPS Features

The current data has 41 subjects, consisting of 3317 total days, and 15 GPS features. To conserve battery life, a subject’s GPS coordinates were tracked for in a 2-min-on and 18-min-off cycle everyday using their own mobile device via the Beiwe platform. We used Ian Barnett’s algorithm to impute the missing data during the off cycles. Assuming no more data was missing due to various factors, one would generate 144 mins of GPS data per day, or 10% of total minutes in a day.

As you can see from the figure below, we can nicely reconstruct an individuals’ mobility trajectory from these data. Figure 1: A weekly view of a subject’s mobility pattern

For an even more intuitive view of the GPS data, take a look at the video here:

From these GPS traces, we extracted daily mobility features, such as max home distance, circadian routine, probability of pauses, as described in Barnett et al., Biostatistics (2020) and defined mathematically in its supplementary material.

To get a flavor of what these features look like, the first six days of GPS data for a subject are attached below.

gps_df = read.table(gps_df_path,header = T, dec = ",", )[,c(1,2,97:111)]

#Define the column types
gps_df$Date = as.Date(gps_df$Date)
gps_df[,3:dim(gps_df)[2]] = apply(gps_df[,3:dim(gps_df)[2]], 2, function(x) as.numeric(x))
head(gps_df)

NA

3. Exclude Data

The first step before further analysis is to exclude data points (days) that had excessive amount of data missing. Here is a historgram of the minutes missing for all 3317 total days.

p=minMiss_histplot(gps_df,200, "All Data")
ggplotly(p)

Figure 2: Minutes missing for all collected days. There are 3317 total days. Dashlines indicate the corresponding percentile.

3a. remove first and last days

First, we will remove the first and last days from each subject, because the application was installed during mid-day at the beginning of the study and uninstalled mid-day at the end of the study.

# loop through each subj to remove 1st and last days of gps data
gps_df_clean = data.frame() #initiate a df
for (subj in unique(gps_df$IID)){ #loop through each subj
  gps_df_subj <- subset(gps_df, IID == subj) #get gps_df per subject
  gps_df_subj <- gps_df_subj[2:(dim(gps_df_subj)[1]-1),] #remove the 1st and last days
  gps_df_clean <- rbind(gps_df_clean,gps_df_subj) #combine all subjs
}

This step removed 82 days. Now, dataset has 41 subjects, consisting of 3235 total days, and 15 GPS features.

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls

Figure 3: Minutes missing after removing the first and last days of each subject. There are 3235 total days. Dashlines indicate the corresponding percentile.

3b. remove days at the sensitivity threshold

Next, we are removing the days with excessive data missingness. This is of course an arbitrary step. Therefore we will conduct a sensitivity analysis of the missingness threshold we set by performing the analysis across multiple different thresholds. For now, the current sensitivity cutoff is set at 1440, which amounts to only excluding those with 100% of GPS missing.

sensitivity_cutoff = 1440 # this controls the cutoff threshold
gps_df_clean2 = subset(gps_df_clean, MinsMissing < sensitivity_cutoff)

This step removed 79 days. Now, dataset has 41 subjects, consisting of 3156 total days, and 16 GPS features.

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls

Figure 4: Minutes missing after removing days with all data missing. There are 3156 total days. Dashlines indicate the corresponding percentile.

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls

Figure 5: A Zoom-in plot of minutes missing distribution after removing days with all data missing. There are 3156 total days. Dashlines indicate the corresponding percentile.

4. Random Data Partitions

To operationalize individual footprinting prediction, we randomly split each individual’s available GPS data to two half partitions for 100 times. We used the features from the first half to determine if we can identify the same individual using their data in the second unseen half. Again, to avoid any (un)lucky random splits, we repeated data partition 100 times.

set.seed(510)
subj_seq = list()
part_times = 100
for (subj in unique(gps_df_clean2$IID)){
  subj_data = subset(gps_df_clean2, IID==subj)
  subj_seq[[subj]] <-createDataPartition(subj_data$IID,times = part_times, p =0.5)
}

5. Build Correlation Matrix

Here, similar to Finn et al. Nat Neuro (2015) and Kaufmann et al. Nat Neuro (2015), we built a Pearson correlation matrix among 16 available GPS features. Below is an illustrative sample of the variables and their correlations.

# an example of subj 1, and first half

Warning messages:
1: In readChar(file, size, TRUE) : truncating string with embedded nuls
2: In readChar(file, size, TRUE) : truncating string with embedded nuls
3: In readChar(file, size, TRUE) : truncating string with embedded nuls
4: In readChar(file, size, TRUE) : truncating string with embedded nuls

example_data = gps_df_clean2[subj_seq$`14w5qlo8`$Resample001,3:17]
gps_cor = rquery.cormat(example_data, type = "full")

Figure 6: Correlation matrix of GPS features for a subject.

We then calculated the feature matrix for all 41 subjects in the dataset, seperately for each half data partition, and for each random split.

make_feature_matrix = function(gps_df) {
  subj_mat_1 = list()
  subj_mat_2 = list()
  for (subj in unique(gps_df$IID)){
    subj_data = subset(gps_df_clean2, IID==subj)
    subj_mat_1[[subj]] = lapply(subj_seq[[subj]], function(list) rquery.cormat(subj_data[list,3:17], type = "flatten", graph = F)$r)
    subj_mat_2[[subj]] = lapply(subj_seq[[subj]], function(list) rquery.cormat(subj_data[-list,3:17], type = "flatten", graph = F)$r)
  }
  return(list(subj_mat_1 = subj_mat_1, subj_mat_2 = subj_mat_2 ))
}

gps_clean2_feature = make_feature_matrix(gps_df_clean2)

6. Match Target

Next we tried to match each target, defined by one subject’s feature matrix in the first half, to a feature matrix in the second half in the same random split. If the same individual’s feature matrix in the second half had the highest correlation among all subjects, then we assigned the match result 1, indicating a succussful match, otherwise 0, indicating an unsccussful match.

# creating a "database" against which target is to be matched with
subj_mat_1 = gps_clean2_feature$subj_mat_1
subj_mat_2 = gps_clean2_feature$subj_mat_2

calc_match_cor = function(subj_mat_1,subj_mat_2) {
  database = list() 
  for (time in 1:part_times) {
      database[[time]] = lapply(subj_mat_2, function(subjmat) subjmat[[time]]$cor)
  }
  
    # match target to database
    match_cor = list()
    for (subj1 in names(subj_mat_1)){ #loop through each subj
      # create a list of target across partitions
      target_list = lapply(subj_mat_1[[subj1]], function(part) part$cor)
      # create a match list
      for (time in 1:part_times){
        target_subj_time = target_list[[time]] #loop through each partition
        for (subj2 in names(subj_mat_2)){ #loop everyone in 2nd half
          data_subj_time = subj_mat_2[[subj2]][[time]]$cor
          match_cor[[subj1]][[as.character(time)]][[subj2]] = cor(target_subj_time,data_subj_time,use = "na.or.complete")
        }
      }
    }
  return(match_cor)
}

match_cor = calc_match_cor(subj_mat_1,subj_mat_2)

By tallying up all the successful and unsuccessful matchs in each of the 100 random splits, we calculated a distribution of match accuracy.

calc_acc_time=function(match_cor,part_times) {
  acc_time = array()
  for (time in 1:part_times){
    acc_time[time] = 0
    for (subj in names(subj_mat_1)){
      max_position = which.max(unlist(match_cor[[subj]][[as.character(time)]]))
      predicted_subj = names(subj_mat_1)[max_position]
      if (predicted_subj == subj) {
        acc_time[time] = acc_time[time] + 1
      }
    }
  }
  acc_time = acc_time/length(names(subj_mat_1))
  return(acc_time)
}
acc_time = calc_acc_time(match_cor,part_times)

Over the 100 random splits, the mean match accuracy was 55.02%, with a standard deviation of 5.9%, high of 70.73%, and low of 36.59%, (95%CI: 0.54-0.56).

p_time_cor = hist_chx(acc_time, bins = 15, title = paste("Correlated Features: \n Average Accuracy Across",part_times,"Data Partitions"), xaxis = "Prediction Accuracy", yaxis = "Count")
ggplotly(p_time_cor)

Figure 7: Match accuracy distributon across 100 data splits. Mean match accuracy was 0.55 (95%CI: 0.54-0.56)

calc_acc_subj = function(match_cor, part_times){
  acc_subj = array()
  for (subj in names(subj_mat_1)){
    acc_subj[subj] = 0
    for (time in 1:part_times){
      max_position = which.max(unlist(match_cor[[subj]][[as.character(time)]]))
      predicted_subj = names(subj_mat_1)[max_position]
      if (predicted_subj == subj) {
        acc_subj[subj] = acc_subj[subj] + 1
      }
    }
  }
  acc_subj = acc_subj/part_times
  acc_subj = acc_subj[-1]
  return(acc_subj)
}
acc_subj = calc_acc_subj(match_cor, part_times)

In addition to the match accuracy across subjects, we also calculated the match accuracy for each individual. Across the 41 subjects, the mean match accuracy was 55.02%, with a standard deviation of 30.19%, high of 98%, and low of 4%.

subj_scatter  = function(gps_df, acc_subj_vector, method){
  subj_df <- data.frame(x=unique(gps_df$IID))
  subj_df$y = acc_subj_vector
  subj_df = subj_df[order(subj_df$y),]
  subj_acc_plot = ggplot(subj_df) + 
  geom_point(aes(x = reorder(x, y), y = y)) + 
   theme_cowplot() + 
    labs(title = paste(method,"Features \n Subject Level Accuracy Across",part_times,"Data Partitions"), 
       x = "Subjects", y = "Prediction Accuracy") +
   theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
  return(subj_acc_plot)
}

ggplotly(subj_scatter(gps_df_clean2,acc_subj, "Correlated"))

Figure 8: Match accuracy distributon across 41 subjects. Mean match accuracy was 55.02%, with a standard deviation of 30.19%, high of 98%, and low of 4%. This shows dramatic differences in individual footprinting distinctiveness.

7. Permutation Test

To assess the statistical significance of the results above, we conducted a non-parametric permutation test. To do this, we randomly scrambled pair-wise subject-to-day linkage. We repeated this process 1000 times, and followed the same procedure as above to obtain a null distribution of average match accuracy across sample, and for each individual.

gps_df_perm = gps_df_clean2
perm_time = 1000
perm_acc_time = list()
perm_acc_subj = list()
for (i in 1:perm_time) {
  perm_part_times = 1
  print(paste("processing ...", i,"..."))
  gps_df_perm$IID = sample(gps_df_perm$IID)
  perm_subj_seq = make_subj_seq(gps_df_perm,part_times = perm_part_times)
  perm_gps = make_feature_matrix(gps_df_perm,perm_subj_seq)
  perm_mat_1 = perm_gps$subj_mat_1
  perm_mat_2 = perm_gps$subj_mat_2
  perm_match_cor = calc_match_cor(perm_mat_1,perm_mat_2)
  perm_acc_time[[i]] = calc_acc_time(perm_match_cor,part_times = perm_part_times)
  perm_acc_subj[[i]] = calc_acc_subj(perm_match_cor,part_times = perm_part_times)
}

Over the 1000 permutations, the mean match accuracy was 2.31%, with a standard deviation of 2.31%, high of 400%, and low of 2.3146341%.(95%CI: 0.02-0.02)

perm_acc_time_all = unlist(perm_acc_time)
q_time_cor = hist_chx(perm_acc_time_all, bins = 8, title = paste("Average Accuracy Across",length(perm_acc_time_all),"Permutations"), xaxis = "Prediction Accuracy", yaxis = "Count")
ggplotly(q_time_cor)

Figure 9: Match accuracy distributon across 1000 permutations. Mean match accuracy in permutation was 0.02 (95%CI: 0.02-0.02). This null distribution does not overlapp with the distributon using real data (Fig.7) at all.

We also calculated null distribution of match accuracy for each subject. Over the 1000 permutations, the mean match accuracy was 2.31%, with a standard deviation of 0.74%, high of 3.7%, and low of 0.6%.

perm_subj = list()
for (subj in subj_acc_plot$data$x) {
  perm_subj$val[[subj]] = sapply(perm_acc_subj, function(perm) perm[which(names(perm) == subj)])
  perm_subj$hist[[subj]] = hist_chx(perm_subj[[subj]]$val, bins = 8, title = paste(subj,": Accuracy Across \n",perm_time,"Permutations"), xaxis = "Prediction Accuracy", yaxis = "Count")
  perm_subj$acc[[subj]] = sum(perm_subj$val[[subj]])/perm_time
}

subj_df <- data.frame(x=unique(gps_df_clean2$IID))
subj_df$y = acc_subj
subj_df$y_perm = unlist(perm_subj$acc)
subj_df = subj_df[order(subj_df$y),]
p_subj_cor_perm = ggplot(subj_df, aes(x = reorder(x, y), y = value)) + 
  geom_point(aes(y = y, col = "subject data")) + 
  geom_point(aes(y = y_perm, col = "permutation")) +
  theme_cowplot() + 
  labs(title = paste("Cor Features \n Subject Level Accuracy Across",part_times,"Data Partitions"), 
       x = "Subjects", y = "Prediction Accuracy") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
ggplotly(p_subj_cor_perm)

Figure 8: Match accuracy null distributon across 41 subjects. Mean match accuracy in permutation was 2.31%, with a standard deviation of 0.74%, high of 3.7%, and low of 0.6%. While individuals exhibited marked differences in footprinting distinctivenss, the subject with the lowest prediction accuracy was still statistically significant against the permutation test (i.e. for subject bkq7uyp6: data: 4 % vs. permutation: 0.6%).

8. Mean Features

We calculated the mean of each GPS feature as a metric of feature variability. Similar to above, we calculated mean separately for each data partition and for each individual.

mean_gps = function(gps_df){
  mean_list = sapply(3:17,function(i) mean(gps_df[,i],na.rm=T))
  names(mean_list) = colnames(gps_df)[3:17]
  return(mean_list)
}

make_feature_mean = function(gps_df, subj_seq) {
  subj_mat_1 = list()
  subj_mat_2 = list()
  for (subj in unique(gps_df$IID)){
    subj_data = subset(gps_df_clean2, IID==subj)
    subj_mat_1[[subj]] = lapply(subj_seq[[subj]], function(list) mean_gps(subj_data[list,]))
    subj_mat_2[[subj]] = lapply(subj_seq[[subj]], function(list) mean_gps(subj_data[-list,]))
  }
  return(list(subj_mat_1 = subj_mat_1, subj_mat_2 = subj_mat_2 ))
}

gps_clean2_mean_feat = make_feature_mean(gps_df_clean2, subj_seq)

Figure 9: Feature mean across all subjects. The histogram are for data in the first random half split.

9. Variability Features

We calculated the root mean square of successive differences or RMSSD of each GPS feature as a metric of feature variability. Similar to above, we calculated RMSSD separately for each data partition and for each individual.

rmssd_gps = function(gps_df){
  rmssd_list = sapply(3:17,function(i) rmssd_id(gps_df[,i], gps_df$IID,long=F))
  names(rmssd_list) = colnames(gps_df)[3:17]
  return(rmssd_list)
}

make_feature_rmssd = function(gps_df, subj_seq) {
  subj_mat_1 = list()
  subj_mat_2 = list()
  for (subj in unique(gps_df$IID)){
    subj_data = subset(gps_df_clean2, IID==subj)
    subj_mat_1[[subj]] = lapply(subj_seq[[subj]], function(list) rmssd_gps(subj_data[list,]))
    subj_mat_2[[subj]] = lapply(subj_seq[[subj]], function(list) rmssd_gps(subj_data[-list,]))
  }
  return(list(subj_mat_1 = subj_mat_1, subj_mat_2 = subj_mat_2 ))
}

gps_clean2_rmssd_feat = make_feature_rmssd(gps_df_clean2, subj_seq)

Figure 10: Feature variability for all subjects. Variability is measured by root mean square of successive differences (RMSSD). The histograms are for data in the first random half split.

10. Predict with New Features

match_mean = calc_match_vector(gps_clean2_mean_feat$subj_mat_1,gps_clean2_mean_feat$subj_mat_2)
acc_time_mean = calc_acc_time(match_mean,part_times)
acc_subj_mean = calc_acc_subj(match_mean, part_times)

match_rmssd = calc_match_vector(gps_clean2_rmssd_feat$subj_mat_1,gps_clean2_rmssd_feat$subj_mat_2)
acc_time_rmssd = calc_acc_time(match_rmssd,part_times)
acc_subj_rmssd = calc_acc_subj(match_rmssd, part_times)

11. Permutation Tests with New Features

11a. Average Accuracy

perm_vector = function(gps_df, perm_time, method){
  gps_df_perm = gps_df
  perm_acc_time = list()
  perm_acc_subj = list()
  for (i in 1:perm_time) {
    perm_part_times = 1
    print(paste("processing ...", i,"..."))
    gps_df_perm$IID = sample(gps_df_perm$IID)
    perm_subj_seq = make_subj_seq(gps_df_perm,part_times = perm_part_times)
    if (method = "rmssd") {

Error: unexpected '=' in:
"    perm_subj_seq = make_subj_seq(gps_df_perm,part_times = perm_part_times)
    if (method ="

#perm_mean = perm_vector(gps_df_clean2, 1000, "mean")


p_time_mean = hist_chx(acc_time_mean, bins = 8, title = paste("Mean Features: \n Average Accuracy Across",length(acc_time_mean),"Data Partitions"), xaxis = "Prediction Accuracy", yaxis = "Count")

q_time_mean =  hist_chx(unlist(perm_mean$perm_acc_time), bins = 8, title = paste("Mean Features: \n Average Accuracy Across",length(perm_mean$perm_acc_time),"Permutations"), xaxis = "Prediction Accuracy", yaxis = "Count")

p_time_mean + q_time_mean

#perm_rmssd = perm_vector(gps_df_clean2, 1000, "rmssd")

p_time_rmssd = hist_chx(acc_time_rmssd, bins = 8, title = paste("RMSSD Features: \n Average Accuracy Across",length(acc_time_rmssd),"Data Partitions"), xaxis = "Prediction Accuracy", yaxis = "Count")

q_time_rmssd = hist_chx(unlist(perm_rmssd$perm_acc_time), bins = 8, title = paste("RMSSD Features: Average Accuracy Across",length(perm_mean$perm_acc_time),"Permutations"), xaxis = "Prediction Accuracy", yaxis = "Count")

p_time_rmssd + q_time_rmssd

combine_perm_subj = function(acc_subj_plot, perm_acc_subj) {
  perm_subj = list()
  for (subj in acc_subj_plot$data$x) {
    perm_subj$val[[subj]] = sapply(perm_acc_subj, function(perm) perm[which(names(perm) == subj)])
    perm_subj$hist[[subj]] = hist_chx(perm_subj[[subj]]$val, bins = 8, title = paste(subj,": Accuracy Across \n",perm_time,"Permutations"), xaxis = "Prediction Accuracy", yaxis = "Count")
    perm_subj$acc[[subj]] = sum(perm_subj$val[[subj]])/perm_time
  }
  return(perm_subj)
}


subj_scatter_perm = function(gps_df,acc_subj,perm_subj, method){
    subj_df <- data.frame(x=unique(gps_df$IID))
  subj_df$y = acc_subj
  subj_df$y_perm = unlist(perm_subj$acc)
  subj_df = subj_df[order(subj_df$y),]
  p = ggplot(subj_df, aes(x = reorder(x, y), y = value)) + 
    geom_point(aes(y = y, col = "subject data")) + 
    geom_point(aes(y = y_perm, col = "permutation")) +
    theme_cowplot() + 
    labs(title = paste(method,"Features \n Subject Level Accuracy Across",part_times,"Data Partitions"), 
         x = "Subjects", y = "Prediction Accuracy") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8), plot.title = element_text(hjust = 0.5))
  return(p)
}

p_subj_mean = subj_scatter(gps_df_clean2, acc_subj_mean, "Mean")
perm_subj_mean = combine_perm_subj(p_subj_mean,perm_mean$perm_acc_subj)
p_subj_mean_perm = subj_scatter_perm(gps_df_clean2, acc_subj_mean, perm_subj_mean, "Mean")
ggplotly(p_subj_mean_perm)

p_subj_rmssd = subj_scatter(gps_df_clean2, acc_subj_rmssd, "RMSSD")
perm_subj_rmssd = combine_perm_subj(p_subj_rmssd,perm_rmssd$perm_acc_subj)
p_subj_rmssd_perm = subj_scatter_perm(gps_df_clean2, acc_subj_rmssd, perm_subj_rmssd, "RMSSD")
ggplotly(p_subj_rmssd_perm)

12. Compare Features

p_time_cor + p_time_mean + p_time_rmssd

p_subj_cor_perm / p_subj_mean_perm / p_subj_rmssd_perm

13. Combine Features

14. Misc Info

#sessionInfo()

 

---

A work by Cedric Huchuan Xia

hxia@upenn.edu